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The Eddy Experiment: GNSS-R speculometry for directional 
sea-roughness retrieval from low altitude aircraft 
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We report on the retrieval of directional sea surface rough- 
ness, in terms of its full directional mean square slope (in- 
cluding direction and isotropy), from Global Navigation 
Satellite System Reflections (GNSS-R) Delay-Doppler-Map 
(DDM) data collected during an experimental flight at 1 km 
altitude. This study emphasizes the utilization of the entire 
DDM to more precisely infer ocean roughness directional 
parameters. In particular, we argue that the DDM exhibits 
the impact of both roughness and scatterer velocity. Ob- 
tained estimates are analyzed and compared to co-located 
Jason-1 measurements, ECMWF numerical weather model 
outputs and optical data. 



1. Introduction 

Several GNSS constellations and augmentation systems 
are presently operational or under development, including 
the pioneering US Global Positioning System (GPS) and the 
forthcoming European system, Galileo. These all-weather, 
long-term, stable and precise L-band signals can be used for 
bistatic remote sensing of the ocean surface and beyond, an 
emerging concept known as GNSS-R. 
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Figure 1. Example of GPS-R Delay-Doppler Map. 



Among several applications, two have been emphasized 
by the community: sea-surface altimetry (see Ruffini et al. 
[2004] and references therein) and sea-surface "speculom- 
etry" (a term discussed below), related to the statistical 
properties of sea surface gravity wave slopes. Although this 
paper addresses the latter, we note the intrinsic capability of 
GNSS-R for providing long term co-located measurements 
of both surface roughness and sea level with high spatial 
and temporal resolution. As recently demonstrated with 
scatterometer measurements (Chelton et al. [2004]), such a 
capability would help to better quantify the relationship of 
velocities in the upper ocean (driven by wind stress forcing) 
with surface height dynamics. 

Thanks to its passive multistatic character, GNSS-R 
clearly holds the potential to provide an unprecedented 
spatio-temporal sampling of the ocean surface. The ex- 
pected high spatial and temporal measurements can cer- 
tainly serve operational applications. For instance, and to 
follow successful scatterometer measurements, GNSS-R can 
complement ocean winds and wave models. Being rain im- 
mune, such new data could help quantify atmosphere-ocean 
coupling, including momentum and energy fluxes under ex- 
treme conditions for hurricane modeling. Many other sci- 
entific applications can be cited, such as sea surface break- 
ing/whitecapping and gas exchange global characterization, 
important ingredients for understanding the ocean's biogeo- 
chemical response to, and its influence on, climate change. 
Indeed, a very promising approach to quantify CO2 flux is 
to better assess the surface fractional area, which is read- 
ily measurable from surface slope measurements ( Watson et 
al. [1999]). In that case, GNSS-R can provide a more di- 
rect measurement to extract gravity surface slope statistical 
properties and to quantify the role of the ocean in taking up 
increases of CO 2. 

In addition, L-band sea surface roughness data can be 
used to support L-band radiometric missions, such as SMOS 
(Soil Moisture and Ocean Salinity) and AQUARIUS, to both 
quantify and efficiently separate roughness and salinity con- 
tributions to L-band radiometric brightness measurements. 

Inferring sea roughness from GNSS-R data requires (i) 
a parametric description of the sea surface, (ii) an electro- 
magnetic and instrument model for sea-surface scattering at 
L-band and (iii) the choice of a GNSS-R data product to be 
inverted. There is quite an agreement on the two first as- 
pects in the literature. It has been recognized that the scat- 
tering of GNSS signals can be approximated by an effective 
Geometric Optics model, where the fundamental physical 
process is the scattering from facet-like surface elements. 
This is the reason for the use of the term "speculometry" 
here, which stems from the Latin word for mirror, speculo: 
the detected GNSS-R return is dominated by the statistics 
of facet slopes and their curvatures at scales larger than the 
electromagnetic wavelength (A). 

Under a Gaussian assumption, three parameters fully de- 
fine the detected L-band sea surface slope probability dis- 
tribution (PDF) . These parameters are encapsulated by the 
directional mean square slope, DMSSa, a symmetric tensor 
which results from the integration of the ocean energy spec- 
trum at wavelengths larger than A, and which characterizes 
the ellipsoidal shape of the slope PDF. 
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Figure 2. DMSSa estimation with the DDM least- 
square inversion approach, along the descending (North 
to South) track— MSS top, SPA middle and SPI bottom. 
The total MSS (Ku-Band) measured by Jason-1 and the 
ECMWF wind direction are also shown for comparison. 



It is important to note that DMSSa has rarely been em- 
phasized as the geophysical parameter of interest in the liter- 
ature. Instead, most authors link sea roughness to the near 
surface wind vector, which is thought to be more useful for 
oceanographic and meteorological users. Unfortunately, this 
is somewhat misleading, as the relationship between surface 
wind and roughness is not one-to-one and requires an addi- 
tional modeling layer. The connection between DMSSa and 
wind is affected by other factors (e.g., swell, fetch and degree 
of maturity), as is well known in the altimeter community 
(see Gourrion et al. [2002]). 



Moreover, the product traditionally used for inversion in 
GNSS-R speculometry is the ID delay waveform of the re- 
flected signal amplitude, from which the wind speed is in- 
ferred assuming an isotropic slope PDF (see, e.g., Garrison 
et al. [1998, 2002] Komjathy et al. [2000], or Cardellach 
et al. [2003]). Attempts have also been made to estimate 
the wind direction by fixing the PDF isotropy to some the- 
oretical value (around 0.7) and using at least two satellites 
reflections with different azimuths (see, e.g., Zuffada et al. 
[2000], Armatys et al. [2000] and Garrison et al. [2003]). 

Here, we will work with a product of higher information 
content: the 2D DDM of the reflected signal amplitude. As 
proposed in Ruffini et al. [2000], the provision of an extra 
dimension opens the possibility of performing a robust es- 
timation of all the DMSSa parameters through the direct 
fitting/estimation of the entire DDM. In Elfouhaily et al. 
[2002], a first order approximation for inversion was pro- 
posed. 

The full DDM-inversion technique proposed in Ruffini et 
al. [2000] is used for the first time here to analyze GNSS- 
R data collected during the Eddy Experiment. This cam- 
paign and the altimetric data analysis is reported elsewhere 
(Ruffini et al. [2004] and Soulat [2003]), and is only briefly 
described in Section 2. The retrieval methodology relies on 
a least-squares fit of the speculometric model, as discussed 
in the third section. In the fourth section, results are com- 
pared to ancillary data (Jason-1 radar altimeter, ECMWF 
numerical weather model, optical data). Finally, another 
important outcome of the exhaustive exploitation of the in- 
formation contained in the DDM product is related to the 
expected mean sea surface motion. We evidence that a small 
part of the Doppler spread can be attributed to the mean 
"scatterer velocity" , i.e., the rapid motion of the A (or larger) 
sized sea-surface facets that contribute the most to the de- 
tected signals. 

2. Data collection and pre-processing 

The data set (i.e., the recorded direct and reflected GPS 
signals together with the aircraft kinematic data) was gath- 
ered during an airborne campaign carried out in Septem- 
ber 2002. The aircraft overflew the Mediterranean Sea, off 
the coast of Catalonia (Spain), northwards from the city of 
Barcelona for about 150 km at 1000 m altitude and 45-75 
m/s speed. The area is crossed by the ground track #187 
of the Jason-1 radar altimeter, which the aircraft overflew 
during the satellite overpass for precise comparison. The 
track was overflown twice: the first time during the ascend- 
ing pass (from South to North) at low speed (45-60 m/s) 
and the second during the descending pass (from North to 
South) at a faster speed (65-75 m/s) due to wind. The time 
shift between the two passes over a same point on the track 
ranged from 45 min to 2h 15 min. During the ascending 
pass, PRNs 08, 10 and 24 were visible with elevations span- 
ning 30° to 85° while PRNs 08, 10 and 29 were visible during 
the descending track with elevations between 40° and 75°. 
The configuration of this test flight was not optimized for 
speculometry: from such low altitude, the sea-surface reflec- 
tive area is essentially limited by the PRN C/A code, and 
the glistening zone is coarsely Delay-Doppler mapped. 

The raw GPS signals were acquired with a modified Tur- 
boRogue receiver, sampled at 20.456 MHz and pre-processed 
with a dedicated software composed of two sub-units fed 
with the direct and reflected signals. Correlations were 
computed at 81 delay lags while the Doppler dimension 
spanned -200 to 200 Hz with a step of 20 Hz. The co- 
herent/incoherent integration times were respectively set to 
20 ms and 10 s, meaning that the averaged DDM were pro- 
duced at the rate of 0.1 Hz after summation of 500 incoher- 
ent looks (see Figure 1 for a sample DDM). 



GERMAIN ET AL.: GNSS-R SPECULOMETRY FROM LOW AIRCRAFT 



X- 3 



3. Speculometric model and DDM inversion 

Following Zavorotny et al. [2000], the link between the 
DDM mean power at delay-Doppler P(t, /) and the effective 
L-band sea-surface slope PDF, P(s x ,s y ), is given by 



P(r, /) = / dxdy 



Cjr r 



J?2 p2 



4-f 

ql 



-qx 
qz 



X [Tm(x,y)-T c -T, fm(x,y)- f c - f], (1) 

where G r is the receiver antenna pattern, R t and R r the 
distances from generic point on sea-surface to transmitter 
and receiver, (q x ,q y ,qz) the scattering vector, \ the Wood- 
ward Ambiguity Function (WAF), T m (x, y) and f m {x, y) the 
delay-Doppler coordinates on the sea-surface and (r c , / c ) the 
delay/Doppler offset of the geometric specular-point with 
respect to the direct signal (the DDM "center"). Account- 
ing for the receiver mean thermal noise Pat and including 
a scaling parameter a, the mean amplitude of the DDM 
can be written as A(r, f) — Wa P(r, /) + Pjv . As discussed 
above, V is described by the DMSSa parameter set, which 
defines an elliptic quadratic form in the 2D space of facet 
slopes. Mean-square slopes along major and minor principal 
axes are often called MSS up-wind (mss u ) and MSS cross- 
wind (msSc) respectively. In the following, we will refer 
to the Total MSS (MSStot~2^/mss u .mss c , proportional to 
the ellipse area and directly related to nadir a"), the Slope 
PDF azimuth (SPA, the direction of semi-major axis with 
respect to North) and the Slope PDF Isotropy (SPI, equal 
to mssc/mss u ). 

The inversion was performed through minimization of the 
mean square difference between model and data DDMs. Nu- 
merical optimization was carried out by a steepest-slope- 
descent algorithm with a Levenberg-Marquardt type adjust- 
ment. The main difficulty stemmed from the presence of 
several nuisance parameters in the forward model (mainly 
r c and f c but also a). The DDM centers were affected by the 
aircraft trajectory (altitude and vertical velocity) to first or- 
der but also by geophysical parameters (such as sea level). 
They needed to be accurately known in order to estimate 
DMSSa- For this reason, the DMSSa and nuisance param- 
eters were jointly estimated in an iterative manner. 

4. Results and analysis 



apparent discrepancy in the southern part of the track can 
be explained as a degenerate solution of the estimation prob- 
lem. Indeed, the inversion of DDM for the SPA is degenerate 
in at least two cases: when the transmitter is at zenith or 
when the receiver moves towards the transmitter. In these 
scenarios, the Delay-Doppler lines mapping the glistening 
zone are fully symmetric around the receiver direction and 
it becomes impossible to distinguish a slope PDF from its 
mirror image about the receiver direction. This effect is 
clearly observed here where the two found SPA (-20° and 
80°) are indeed symmetric around the aircraft heading direc- 
tion (30°). In this part of the track, the azimuth of PRN-08 
is about 50°, almost aligned with the aircraft heading di- 
rection. In the northern part of the track, the estimated 
SPA matches very well with the wind direction provided by 
ECMWF. In the southern part, the mismatch reaches up 
to 50°. However, we underline that surface wind is not the 
only element driving SPA and that swell is likely to have 
contributed. 

Finally, the bottom plot shows the SPI variations along 
the track. The reflected signals are strongly directional. The 
wind-driven spectrum of Elfouhaily et al. [1997] for a ma- 
ture sea predicts a SPI value around 0.65, largely insensitive 
to wind speed. The apparent significant departure from this 
reference value is a probable signature of an under-developed 
sea with the presence of swell. However, the relatively poor 
consistency among PRN remains an issue to clarify: further 
work is needed to validate the accuracy of these SPI esti- 
mates and to better understand the potential information 
(and possible applicability constraints) of this product. 

As a second and new outcome of this analysis, we now 
discuss the signature of "scatterer velocity" in the data, i.e., 
the signature of fastly moving sea-surface facets with size 
(curvature) larger than ~20 cm. Such a signature can be 
detected when comparing the total MSS along the ascend- 
ing and descending tracks, as estimated by the least-squares 
approach: a drastic discrepancy (up to 33%) was observed 
for two passes shifted by less than hour over the same track 
point. 

Multipath effects could conceivably lead to a Doppler 
width modulation. However, some azimuth dependence 
should have been observed and was not. Another possible 
cause could be a changing aircraft attitude between ascend- 
ing and descending tracks, but the aircraft roll and pitch 



The values of DMSSa estimated along the descending 
track of the flight are shown on Figure 2. The top plot 
illustrates the variations of Total MSS. The inter-PRN con- 
sistency is reasonable in the southern part of the track but 
worsens slightly in the northern part. For comparison, the 
total MSS in Ku-band was derived from the Jason-1 a co- 
located measurements at 1 Hz sampling (7 km) and 20 km 
resolution. The Jason-1 MSS was obtained through the sim- 
ple relationship MSS= k/<j°, k being the effective (empir- 
ical) Fresnel coefficient, here set to 0.45. As expected, we 
observed that the level and dynamic of MSS decreased with 
longer wavelength (from 2 cm in Ku-band to 19 cm in L- 
band). The lower dynamic of L-band MSS impeded any 
clear trend comparison, although the measurements of PRN- 
10 seem in good agreement with Jason-1. The Jason-1 wind 
speed, derived from both the Ku-band cr° and the signifi- 
cant wave height (of about 2 m), ranged from 9 to 13 m/s 
along the track. Translating this wind speed into L-band 
MSS through the spectrum of Elfouhaily et al. [1997] yields 
values between 0.0220 and 0.0255, in-line with GNSS-R re- 
sults. However, we must emphasize that the assumption of 
a wind-driven spectrum was not really warranted during the 
campaign. 

SPA estimation results are presented on the middle plot. 
The inter-PRN consistency is here very satisfying, and the 
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Figure 3. Average scatterer velocity obtained when 
assuming a perfect match of ascending/descending MSS 
and the first order MSS model. It correlates fairly well 
with wind speed and the observed swell from optical data. 
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values were checked to be nominal along both tracks. More- 
over, a change in yaw would slightly impact the antenna 
pattern ground projection but would translate only into a 
Doppler bandwidth cut and never a broadening. We con- 
cluded that the most likely explanation was of geophysical 
origin: the Doppler spectral width had been modulated by 
the overall sea surface dynamics as detected from the mean 
specular facet motion. While the inversion approach as- 
sumed a still surface, the relative velocity between receiver 
and scatterer should be taken into account for proper DDM 
inversion. At high receiver velocities this assumption is fine 
because the scatterer velocity impact will not be significant. 
At low speeds, however, scatterer velocity becomes relevant. 
This analysis is consistent with the fact that the MSS es- 
timated in the ascending track (slower aircraft speed, from 
45 to 60 m/s) showed abnormal high values compared to 
the ones estimated during the ascending track (faster speed, 
from 65 to 75 m/s). 

In order to test this idea in a simple manner, we used 
the results in Elfouhaily et al. [2002], where it was demon- 
strated that a first order relationship exists between the 
moments of the DDM and the full set of DMSSa if the im- 
pact of the bistatic WAF and antenna gain are neglected: 
MSS=A 2 B 2 /(2V r2 sin 2 e), where e is the transmitter eleva- 
tion, V the receiver speed and B the DDM Doppler band- 
width. Assuming that the MSS did not vary significantly be- 
tween the times of the ascending and descending passes and 
considering the relative velocity between sea-surface scat- 
terers and aircraft, i.e., V ±v a for the descending/ascending 
passes respectively, we can roughly solve for the scatterer ve- 
locity v s (for simplicity assumed here parallel to the ground 
track). Applying this scheme to all possible pairs of mea- 
surements and averaging, the plot of Figure 3 results. The 
apparent mean scatterer velocity variations correlate with 
the JASON-1 wind speed variations. Moreover, from the 
dispersion relation an average scatterer speed of 8 m/s cor- 
responds to waves of about 45 m wavelength. The detected 
motion thus seems to be associated with the longer and more 
coherent wave components, consistent with optical obser- 
vations of the swell vector (wavelength and direction, see 
Soulat [2003]), revealing the presence of a northerly (almost 
aligned with the flight track) generated 48 m wave system. 
As hypothesized, the quasi-specular facets contributing the 
most to the bistatic reflected signals may, on average, be 
assumed to almost coherently travel with the dominating 
longer wave peak component. 

5. Conclusion 

We have reported the inversion of GNSS-R signals us- 
ing, for the first time, the full Delay-Doppler Maps for the 
retrieval of the sea-surface directional mean square slope, 
DMSSa: the estimates show good inter-PRN consistency 
(except for the measured anisotropy SPI) and fair agree- 
ment with other sources of data. 

The use of the full DDM further has helped to reveal a 
geophysical signature in GNSS-R associated with the mean 
sea surface scatterer velocity. Under quasi-specular condi- 
tions, sea surface scatterers are mostly associated to small 
slopes corresponding to longer waves with velocities that can 
reach 5-10 m/s, impacting significantly the Doppler band- 
width of slow-moving receivers (e.g., airborne or ground- 
based, Soulat [2004]). The detection of such a geophysical 
signature opens new opportunities for GNSS-R speculome- 
try: to infer either DMSSa or a combination of DMSSa and 
scatterer velocity, depending on the aircraft speed. Further 
investigations will be carried out to take into account the 
correct deformation of the Doppler lines on the surface for 
the search of a scatterer velocity vector. Finally, we em- 
phasize that the flight was not optimized for speculometry: 
higher and faster flights are needed in the future to consol- 
idate the DDM inversion technique and to test new higher 
resolution inversion concepts. 
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